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Abstract 

Generation of thermal convection flow in the liqnid metal battery, a device 
recently proposed as a promising solntion for the problem of the short-term 
energy storage, is analyzed nsing a nnmerical model. It is fonnd that con¬ 
vection cansed by Jonle heating of electrolyte dnring charging or discharging 
is virtnally nnavoidable. It exists in laboratory prototypes larger than a few 
cm in size and shonld become mnch stronger in larger-scale batteries. The 
phenomenon needs fnrther investigation in view of its positive (enhanced 
mixing of reactants) and negative (loss of efficiency and possible disrnption 
of operation dne to the flow-indnced deformation of the electrolyte layer) 
effects. 
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1. Introduction 

The work reported in this paper is motivated by the recent attempts to 
develop a commercially viable liqnid metal battery, a device for large-scale, 
short-term, stationary energy storage [H, . The battery, originally proposed 

in 1960s (see, e.g., |^), is now a snbject of renewed attention as a promising 
solntion of the problem of intermittency of energy snpply from wind and 
solar sonrces j^. 

The energy stored in a liqnid metal battery is the difference between 
the Gibbs free energy of a free light metal (e.g. Na, Li, or Mg) and of the 
same metal in a componnd with a heavy metal (e.g. Bi, Sb, or PbSb). The 
processes of charging or discharging the battery correspond, respectively, to 
electrochemical redaction of the light metal from the componnd or forming 
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Figure 1: (a)^ Model geometry of the battery’s model and the coordinate system. Three 
liquid layers B, E, and A fill a cylindrical cavity. During the charging and discharging pro¬ 
cesses the imposed uniform vertical electrical current of density Jq = Jo^z generates inter¬ 
nal Joule heating and the azimuthal magnetic field Bq = B^eg. (b), Typical distribution of 
non-dimensional temperature in the model battery in the absence of convection-generated 
flow (see section lOI for an explanation). 


the compound. The reactions happen entirely in liquid state at the interfaces 
between the metals and a molten-salt electrolyte, which separates the metals 
from each other, immiscible with them, and is conductive to the ions of the 
light metal. An example of the electrolyte is LiF-LiCl-Lil, which can be used 
in a battery operating at tenmeratures about 450°C with Li as a light metal 
and SbPb as a heavy metal [j]. 

The possible combinations of materials and features of the currently pur¬ 
sued battery designs are discussed, for example, in and j^. Important for 
us is that an operating battery can be viewed, in a simplihed way, as a sys¬ 
tem schematically represented in Fig. [TK, i.e., as a cylindrical box hlled with 
three layers: the bottom layer B containing mixture of the heavy metal and 
the compound, the top layer A containing the light metal, and the electrolyte 
layer E in the middle. The system is stably stratihed, with the densities of 
the materials satisfying ps > Pe > Pa- During the charging and discharging, 
strong (about 1 A/cm^) electric current flows between the top and bottom 
walls serving as current collectors. The sidewalls are electrically insulating. 

The starting point of our work is the recognition that the system can 
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experience hydrodynamic instabilities leading to flows in all three layers. The 
flows would signiflcantly affect the operation of the battery, with implications, 
which can be both positive (faster electrochemical reactions due to enhanced 
mixing in the layer B) and negative (deformation of interfaces leading to 
nonuniform reaction rates and, in extreme cases, short circuit between the 
layers A and B). 

The instabilities can be considered as an essential part of the general scale- 
up problem, i.e. the problem of transition from small laboratory prototypes 
to large, commercially viable devices. It appears inevitable that increased 
size would lead to new instabilities, which would change the hydrodynamics 
of the system and affect its operation. 

The hydrodynamics of a liquid metal battery is largely unexplored, but 
general physical reasoning and preliminary studies suggest existence of sev¬ 
eral distinct instability mechanisms. One is the Tayler instability, the special 
case of the non-axisymmetric pinch-type instability in a column of an electri¬ 
cally conducting fluid with axial electric current. The instability, tradition¬ 
ally considered in the astrophysical context, has been recently analyzed for 
single-metal columns representing the battery in a drastically simplified way 
js], 0,0,0] • It has been found that the instability threshold is predominantly 
determined by one non-dimensional parameter, the Hartmann number 

Ha = Bo(R)R.I^, (1) 

V 

where Bo{R) is the magnitude of the azimuthal component of the magnetic 
held created by the axial electrical current at the outer radius R of the 
column, and a, p, v are the electrical conductivity, density, and kinematic 
viscosity of the metal. In the simplest case of an infinite column with the base 
state having zero flow velocity and uniform vertical current, the instability 
first occurs when Ha exceeds Hacr ~ 22 and has the form of exponential 
growth of perturbations with the azimuthal wavenumber m = 1 and the axial 
wavelength about O.Svri? 0]. Estimates based on the physical properties of 
typical liquid metals show that the instability occurs in moderately sized 
batteries. For example, a column of Mg at 750°C with axial current of 1 
A/cm^ is unstable at R above approximately 25 cm. We should note that, 
as proposed in 0], the Tayler instability can be shifted to higher Ha or even 
completely avoided via modifications of battery’s design that alter the base 
magnetic held Bq. Furthermore, as found in the recent analysis [l^, even 
when the instability occurs in moderately sized batteries, the amplitude of 
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the developing flow is small. Its kinetic energy is insnfficient to overcome the 
gravitational potential of the stably stratified three-layer system and, so, to 

that the Tayler 
^ 1.5 m or even 


1C 


cause rupture of the electrolyte layer. It is suggested in 
instability in a typical battery first becomes dangerous at R 
higher. 

Existence of another instability mechanism is suggested by the analogy 
between the liquid metal batteries and the aluminum reduction cells, the 
devices in which aluminum is produced from oxide by the electrochemical 
Hall-Heroult process. In a simplified description, a reduction cell is a large 
(about 4 by 12 m horizontally), shallow rectangular bath filled by a layer 
of molten aluminum at the bottom and a layer of molten salt electrolyte, 
with aluminum oxide dissolved in it, at the top. It is known in the indus¬ 
try that the cells may experience instability in the form of growing sloshing 
motions of the interface between the two layers. The instability is magneto¬ 
hydrodynamic in nature, related to the large difference between the electrical 
conductivities of electrolyte and aluminum, and caused by the interaction be¬ 
tween the magnetic held generated by external electrical conductors and the 
perturbations of electrical currents within the cell arising due to the local 
deformations of the interface [ll[ 12 
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In a liquid metal battery, the electrical conductivity of electrolyte is about 
four orders of magnitude lower than the conductivities of both metals. Any 
non-uniformity of the thickness of layer E would cause significant redistribu¬ 
tion of electric currents within the battery and, thus, change of Lorentz forces 
acting in the melts. Since no results concerning this phenomenon have been 
published yet, we have to consider the possibility of the additional Lorentz 
forces leading to such a multi-layer magnetohydrodynamic instability as hy¬ 
pothetical. We should also note that the analogy between the battery and 
the aluminum reduction cell is incomplete, primarily because of significant 
differences in geometry. 

The compositional convection can be caused by spatial variations of the 
concentrations of the heavy metal and compound in the bottom layer. While 
likely to be a factor affecting the battery’s operation, this mechanism is im¬ 
possible to analyze at the moment due to complexity and uncertainty of the 
relevant phase diagrams. We should also mention the short-wave instabil¬ 
ities of either Marangoni or electrohydrodynamic nature that may appear 
at the interfaces between the layers and may play a role in the dynamics of 


the battery. Finally, as suggested recently in [8|, flow in a battery can be 
generated by the electrovortex instability near the current collectors in the 
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top and bottom walls. 

Yet another, virtnally nnavoidable instability mechanism is that of ther¬ 
mal convection cansed by non-nniformity of the temperatnre held. There are 
two sonrces of the non-nniformity: heating of the walls applied to maintain 
the necessary operational temperatnre and the internal Jonle heating of the 
poorly condncting electrolyte. The wall heating is reqnired in small-scale 
laboratory prototypes bnt may become nnnecessary in larger cells, where the 
Jonle heating will be snfhcient or even excessive (so the battery will have to be 
cooled) for temperatnre maintenance. To onr knowledge, the only pnblished 
resnlts on convection in liqnid metal batteries are those of the experiments 
(lij l. in which a simplihed system consisting of a single liqnid metal layer 
heated from below is considered. The focns of this work is on the expected 
positive effect of convection as a mechanism of mixing in the bottom layer. 

Our paper presents the hrst study of the convection caused by the internal 
Joule heating in the electrolyte layer. We apply numerical simulations to 
determine the typical size of the battery, at which convection flows become 
inevitable, and to understand the nature and properties of the flow. 


2. Model and approach to analysis 


2.1. Physical model and governing eguations 

Even in its simplified form illustrated in hguredl a liquid metal battery is 
a complex electro-magneto-hydrodynamic system. Its dynamics is a result of 
combined action of many mechanisms, some of which are possibly unknown. 
A reasonable hrst step of the analysis is to explore the known mechanisms 
individually, using idealized systems, in each of which the action of a par¬ 
ticular mechanism is isolated. This was done, for example, for the Tayler 
instability in B 0, H B [13 by neglecting the layered nature of the system 
and heating, i.e. by considering a column of an isothermal liquid metal. 

Here, we analyze the isolated effect of thermal convection. The effects 
of the Tayler instability, compositional convection, and the Lorentz forces 
associated with deformation of interfaces, short-wave interfacial instabilities, 
and electrovortex instabilities are removed by simplifying assumptions. As we 
show later, the approach is partially justihed by the fact that the convection 
already occurs in batteries of very small size, signihcantly smaller than the 
sizes, at which the Tayler and, likely, other instabilities become active. 

The simplifying assumptions made in our model are as follows: 
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1. The sidewalls of the battery are assumed thermally and electrically 
perfectly insulating. 

2. The top and bottom walls, which, in a real battery, serve as cur¬ 
rent collectors, are assumed to be perfectly electrically conducting and, 
therefore, modeled as equipotential surfaces. We also assume that the 
top and bottom boundaries are maintained at a constant temperature. 
This corresponds to an external heating or cooling arranged in such a 
way that, together with the internal Joule heating, it creates a desired 
steady operational temperature within the battery. 

3. The layered nature of the system is only partially retained in our model. 
Specihcally, the electrical and thermal conductivities are assigned dis¬ 
tinct and realistic values in each layer. For the other physical prop¬ 
erties, such as density, viscosity, specihc heat, and thermal expansion 
coefficient, we use the values typical for the electrolyte in all three lay¬ 
ers. This departure from reality allows us to focus the analysis on the 
convection caused by non-uniform internal Joule heating and affected 
by the magnetic held. It also allows us to use the Boussinesq approx¬ 
imation, which would be impossible at hnite density differences be¬ 


4. 


tween the layers [16|, and to apply the effective computational method 


described in section 12.41 

The base state is that of zero how velocity in all three layers and the 
electrical current of constant and uniform density 


Jo = JoG-z, Jq = const (2) 

howing between the top and bottom walls (see hgured^). The current 
generates the constant, uniform, and purely azimuthal magnetic held 

Bo = Bo{r)ee = (3) 

where po is the magnetic permeability of free space. 

5. The interfaces between the layers are modeled as steady-state, horizon¬ 
tal, and impermeable surfaces. The ehects related to the deformation 
of the interfaces are, thus, neglected in our model. This is justihed by 
the fact that the actual battery systems are strongly density-stratihed. 
The density of the light metal is, at least, two times smaller than the 
density of the electrolyte, which, in turn, is several times smaller than 
the density of the compound between the light and heavy metals. A 
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further discussion of this approximation based on the results of our 
computations and supporting its validity is given in section |U 

6. The coupling between the flows in adjacent layers via viscous shear 
stresses, heat transfer, pressure forces, and electromagnetic effects is 
included into the model. The effect of the surface tension at the inter¬ 
faces is assumed to be weak in comparison to these mechanisms and 
neglected. 

7. The fluid is assumed to be Newtonian, incompressible, electrically con¬ 
ducting and, apart from the stepwise changes of thermal and electrical 
conductivities at the interfaces, having constant physical properties. 

8. We assume that the typical timescale of the evolution of convection 
flows is much smaller than the typical time of charging or discharging 
the battery. For this reason, the thicknesses and physical properties of 
the three layers are assumed constant. 

9. The Boussinesq approximation is used to describe the convection effect. 

10. The quasi-static approximation 0 is used to evaluate the electric 

current perturbations induced by the flow velocity and the Lorentz 
forces resulting from the interaction of these currents with the base 
magnetic held B. The approximation is considered valid because the 
magnetic Reynolds number based on the typical velocity and length 
scales of the convection how is expected to be small. The perturbations 
of the magnetic held and, thus, the mechanism leading to the Tayler 
instability are not considered. 

To non-dimensionalize the governing equations and boundary conditions, 
we use the physical properties of electrolyte: density pe, kinematic viscosity 
ue, electrical conductivity aE, thermal conductivity ke, specihc heat C®, 
and thermal expansion coefficient ag. The radius of the battery R is used 
as the typical length scale. In order to derive the temperature and velocity 
scales, we need, hrst, to consider the internal Joule heating by the imposed 
current Jo- In the electrolyte, the volumetric heating rate is 

Qo = —. (4) 

O'E 

This constant will be used as the typical scale for the internal heating rate, 
which is equal to and in, respectively, layers A and B. The 

typical scales of temperature and velocity can then be set to 

^ ^ JlR^ ^ ^ ^ ^oeqATR, (5) 

Ke Ke 
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where g is the gravity acceleration constant. The typical time and pres¬ 
sure scales are RU~^ and PeU"^- For the density of the additional electrical 
currents induced by the flow in the imposed magnetic held Bq, we use the 
Ohm’s law to derive the typical scale aEUB^. The corresponding scale for 
the electric potential is UBqR. The imposed magnetic held Bq is scaled by 
its magnitude at the sidewall Bq{R) = pqJqR/2. 

The imposed temperature of the top and bottom walls is taken as the 
reference temperature and set to zero. The gravity force and the base-state 
Lorentz force —JoB^er resulting from the interaction between the imposed 
current Jq and the magnetic held Bq are removed from consideration by sub¬ 
tracting from the total pressure held the steady-state pressure distributions 
balancing these forces. 

The non-dimensional governing equations are: 

du , _ 1 Ha^ / . „ s ^ 

-h X Bq) +Te^, 


dt 

V ■ u = 0, 
dT 
dt 


u ■ V)w = —Vp 


1 

Re 


Re 


— + U- VT = 


1 


RePr 

j = a(-V0 + UX Bo), 

V ■ (ctV0) = V ■ {au X B 


[V-{nVT) + Q], 


( 6 ) 

( 7 ) 

( 8 ) 
(9) 

( 10 ) 


where u, T, j, cj), and p are the non-dimensional velocity, temperature, den¬ 
sity of electric currents induced by the how, the electric potential associated 
with these currents, and the modihed pressure. The equations are written 
for the entire how domain shown in hgure [H The non-dimensional rate of 
Joule heat generation Q, thermal conductivity n, and electrical conductivity 
a change discontinuously at the interfaces between the layers and are dehned 
as: 



aX z < Ze — He/^,, 


a = 1 1 

at Ze — He/2, < z < ze R He/2, 

(11) 

( cta/cte 

at z > Ze + He/2,, 


( Kb/ke 

aX z < Ze — He/2,, 


K = < 1 

at Ze — He/2, z < ze -\- He/2,, 

(12) 

[ I^a/I^E 

aX z > Ze + He/2,, 



where ze and He stand for the location of the midplane and thickness of 
the electrolyte layer, so z < ze — He/^,, ze — He/^, < z < ze + He/^,, and 
z > Ze + He/^, correspond to the layers B, E, and A, respectively. 




The boundary conditions are 


(13) 

(14) 


T = 

dT 

dr 


0, u = 
= 0, u 


0, 0 = 0 at z = 0, H, 
dd> 

= 0, ^ = 0 atr = 1. 
or 


Additional boundary conditions follow from our assumption of non-deformable, 
impermeable, horizontal interfaces, which implies zero vertical velocity at 
them: 


Uz = 0 at z = ze ^ He/^. (15) 

The non-dimensional parameters of the problem are the Reynolds number 


where 


Re^ — = Gr^l\ 
ve 


^ aEgQpR^ _ aEgJpR^ 

is the Grashof number, the Prandtl number 

J^ePeCe 
Fr = -, 


Ke 


(16) 

( 17 ) 


(18) 


the Hartmann number 


Ha = Bo{R)R 


o'e 

PeJ^e 


1/2 


PqJq p2 f \ 
2 \PeEe ) 


(19) 


and the geometric parameters, namely the aspect ratio H of the battery, and 
the thickness He and location ze of the electrolyte layer. 

As the base state of the system, we consider the always existing mathe¬ 
matical solution with zero velocity and electric potential, constant modihed 
pressure, and the piecewise-parabolic velocity prohle Tq{z) found as an ana¬ 
lytical solution of 

-Q{z), (20) 

0 , ( 21 ) 

where e{z) and Q{z) are dehned by flTT]) and ([12]). As an illustration, the 
prohle of Tq{z) at ze = 0, He = 0.304 is shown in hguredb. 


d r , .dTo 
Tz 

To{0)=To{H) = 
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2.2. Approach to analysis 

The convection instability and the resnlting flows are analyzed via di¬ 
rect nnmerical solntion of the system of eqnations and bonndary conditions 
Each simulation starts with random low-amplitude (~ 10“^) pertur¬ 
bations of velocity and temperature added to the base state. The base state 
is deemed stable or unstable if, after initial adjustment, the energy of the 
perturbations, respectively, decays or grows over a long (at least a hundred 
units) time period. In all the cases of growing perturbations, we have been 
able to identify sufficiently long periods of exponential growth, during which 
the growth rate 

1 dE 1 dE^ 

of the volume-averaged perturbation energy 

E=^ [ \u\‘^dv OT E^ = ^ [ {T- Tofdv (23) 

y .Jv V Jv 

is constant within the second significant digit. 

The growth of perturbations is followed into the stage of a fully developed 
flow characterized by nonlinear saturation. After that, the flow’s evolution is 
computed for not less than 1000 time units. Flow statistics are accumulated 
and time-averaged during this period. 

2.3. Parameter range under investigation 

As we have already discussed, the convection instability is expected to 
appear in the process of scale-up of a battery from small laboratory proto¬ 
types to large commercial devices. The important questions are: at what 
size the convection first appears in a battery of a certain design, and how the 
strength of convection and its mixing and heat-transfer effects vary with the 
size. We choose to address the questions and, so, vary the non-dimensional 
parameters of the problem in the manner that corresponds to variation of 
the size of a battery of a given design. We also explore the effect of the thick¬ 
ness of the electrolyte layer. The physical properties of the melts, imposed 
electric current density, the aspect ratio of the cell, and the location of the 
mid-plane of the electrolyte layer are assumed constant. This means that we 
are left with two independent parameters: the Grashof number Gr and the 
non-dimensional electrolyte thickness He- The other parameters defined in 
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section HH] are either constant (Pr, H, ze, crB/uE, o'a/o'e, i^a/i^e) 

or fnnctions of Gr\ the already mentioned Re = and 


Ha = [5Gr^^\ 


(24) 


where the coefficient 

JI^eI^e^e \ 

2 V ) 


(25) 


(see (HZD and ([I9D). 

Cells of the aspect ratio H = 1 and with the imposed electric current 
Jo = 1 A/cm^ 0,3,3 are considered. The choice of physical properties is 
more difficult. The available data on the properties of the high-temperature 
melts used in the batteries are incomplete and, generally, of low quality. 
Further uncertainty is created by the presence of the ions of the light metal in 
the electrolyte, mixture of the heavy metal and the compound in the bottom 
layer, and the fact that the related phase diagrams are poorly known. For 
these reasons, we select a set of physical properties, which is typical in the 
sense that it represents the properties of the battery materials on the level 
of the orders of magnitude, but does not correspond to any specific battery 
design. 

The well-documented properties of LiCl-KCL at about 450°C are chosen 
for the electrolyte: pe = 1.63 x 10^ kg/m^, aE = 2.93 x K“^, Ce = 
1.21 X 10^ J/kg-K, ue = 0.71 X 10“® m^/s, ke = 0.42 W/m-K, cte = 170 S/m 
(see 0,0). This gives 


Pr = 3.33, 13 = 1.05 x 10“®. 


(26) 


The battery’s radius is related to the Grashof number as 

R = ( Gr ^^^^ ] = 6.60 X 10 “^ [ml. 

V oieqJo J 

The ratios of electrical and thermal conductivities are set at 
^ ^ = 10 ^ ^ ^ = 50 , 

O'E O'E Ke f^E 


(27) 


(28) 


which corresponds to typical high electrical and thermal conductivities of 
liquid metals. 
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2.4- Numerical method 

The evolution of a three-dimensional unsteady flow is calculated in the 
manner of direct numerical simulation using the method introduced in 19 


and extended to flows with cylindrical geometry in |20|, and flows with 


the combined effects of convection and magnetic held in [2l|, l22|, l23|, l2J, l25 


The details of the method can be found in these references, while only the 
main principles and the new features introduced in this study are described 
here. 

The method is based on the second-order hnite-difference discretization 
on a collocated, structured, non-uniform grid. The spatial derivatives are 
evaluated with the use of velocity and current huxes obtained by interpo¬ 
lation to staggered grid points (see (l9|). This makes the scheme nearly 


fully conservative in the sense that, in the non-dissipative limit, it perfectly 
conserves mass, momentum, electric charge, and internal energy, while the 
kinetic energy is conserved with the dissipative error of the third order. 

The time-discretization uses the standard projection method to hnd pres¬ 
sure and satisfy incompressibility, and the second order backward-difference- 
Adams-Bashfort scheme modihed so that the temperature diffusion term is 
treated implicitly. Each time step is accomplished as a sequence of the fol¬ 
lowing substeps: 

(i) Use the fields obtained at the previous time level to solve the 
potential equation and find the electric currents 


V ■ (aV0”) = V ■ (an" x Bo ), 

r = a(-V0" + n"xe,), 


(29) 

(30) 

(31) 


and compute 


F" = -M(n", n") + ^V'n" + X Bo + Te„ 

Re Re 


= -V ■ (T"n") + 




RePr^ 


(32) 

(33) 


where AT(n",n") is the nonlinear term in divergence form. 
(ii) Find the intermediate velocity field u* from: 


3n*-4n" + n" ^ ^ 
2At 


( 34 ) 
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(iii) Solve the Poisson equation for pressure and correct the velocity held 
to satisfy incompressibility: 

VV+‘ = (36) 

u"+‘ = u" - ?AiVp”+‘. (36) 

(iv) Solve the implicit equation for temperature: 

- ^ -= 2G” - G”-' + —V . (KVr”«). (37) 

The impermeability boundary conditions for velocity are satished via the 
pressure correction as discussed below. The no-slip velocity boundary condi¬ 
tions at solid walls are imposed explicitly after 0361) . The boundary conditions 
for potential and temperature at the walls are satished as a part of solution 
of the elliptic problems fl29l) and fl371) . 

For the computations conducted in this work, the structured grid is built 
on the lines of the cylindrical coordinate system (see hgure [I]) and clus¬ 
tered toward the walls and the interfaces between the layers according to 
the coordinate transformation successfully used in our recent work . The 
clustering is a mixture of the Chebyshev-Gauss-Lobatto clustering and the 
uniform grid. It implements the advantages of a nearly-Chebyshev resolution 
of boundary layers, while avoiding the strict time-step limitations caused by 
the smallest grid step near the boundary. 

In the radial direction, the coordinate transformation is 


r = Ar sin(7r?7/2) + (1 — Ar)rj, (38) 

where 0 < r; < 1, ^4^ = 0.96, and the grid is uniform in the r^-coordinate. 

In the vertical direction, a similar transformation is applied separately to 
each layer in such a way that points are clustered towards both the upper 
and lower boundaries. Denoting the global transformed coordinate, in which 
the grid is uniform, as — 1 < .^ < 1 and introducing local coordinates 

6 = = (39) 

Kt - 6 Zt- Zb 

where Zb, ^b and Zt, are the global coordinates of the bottom and top 
boundaries, we dehne the transformation as 

ze = sin«f/2) -|- (1 - A^)^e. (40) 
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The blending coefficient is = 0.96 except for the case of thin electrolyte 
layer He = 0.1, in which the grid is uniform, so = 0. 

The total number of grid points in the vertical direction is divided equally 
among the three layers. 

Two major new features have been introduced into the algorithm for the 
purposes of this work. One of them concerns the velocity boundary condition, 
which now includes impermeability at not just solid walls, but also at the 
interfaces between the layers (see In order to enforce that we solve the 

pressure equation (135|) separately in each of the three layers A, E, and B and 
impose the Neumann condition obtained by projection of fl36D on the normal 
to the boundary 


- = - 1 ! 

On 2At ” 

at the entire boundary of each domain, including the interfaces. The pressure 
fields obtained in this solution are dehned up to additive constants, and can 
be combined into one continuous pressure held by adjusting the constants. 

Another new feature concerns the solution of the elliptic problems (|2^ 
and (jST]) for electric potential and temperature. In the original and new 
versions of the algorithm, this is done using the Fast Fourier Transform in 
the azimuthal angle 6 and the direct solution of the two-dimensional elliptic 
problems for the Fourier coefficients - functions of r and z. Implementation 
of this procedure for fl29l) and fl37|) applied to the entire how domain requires 
modihcations accounting for the presence of variable coefficients a{z) and 
k{z). As discussed in detail in [l9|, the 2D elliptic problems for the Fourier 
coefficients are solved using the high-level routine of the FishPack library 
that discretizes a 2D separable elliptic equation on a uniform structured grid 
using central diherences of the second order and solves the resulting matrix 
equation by the cyclic reduction method. In order to use this algorithm and 
retain the conservation properties of our scheme, we have to express the 2D 
elliptic equations in the transformed coordinates t]-^ in such a way that the 
discretized equations produced by the FishPack routine would correspond 
exactly to the discretization of the original equations on the non-uniform 
grid on the r- 2 ;-plane. 

The procedure for the equations with constant coefficients is described 
in [l^. For the case of variable coefficients, it is modihed as illustrated by 
the following example. We divide the potential equation fl2^ by a{z) and 
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require that the discretization on the non-uniform grid: 


1 5 

a 6z 


a 


Scf) 

6z 


1 


a. 


3 1/2 z —; 


pj-1 


■J-1 


a,- 


Zj+l/2 - Zj-ll2 


(42) 


is reproduced exactly by the central-difference approximation on the uniform 
grid of the same term expressed in the transformed coordinates (we use the 
same notation as in 0 ): 


1 5 
a 6z 


a 


6(p 

6z 




20 j -f 0 j_i 4>j+i ~ ~ 1 




C 2 - 


2Ae 


(43) 


In the formulas above, 6 stands for the central-difference operator, and the 
quantities at the half-integer points are obtained by interpolation, e.g., as 
Zj +\/2 = {zj+i + Zj)/2. Direct comparison of the two formulas leads us to: 


Cl = 

(Zj+l/2{,Zj - Zj-i) + aj-i/2{Zj+i - Zj) 

(44) 

2A^ 

C2 = 

^■^-1/2(^ 1+1 ~ Zj) + aj^i/2{Zj — Zj-i) 

(45) 

Ae 

C 3 = 

%'+i ~ ^'-1 ~ ^3 ~ 

' 2Ae Ae Ae ■ 

(46) 


The other terms of the equations involving the variable coefficients a{z) 
and k{z) are discretized directly on the non-uniform grid in a straightforward 
manner. 

The algorithm is parallelized using the OpenMP parallelization, with a 
typical simulation ran on 16 cpus of a shared-memory workstation. 

A grid sensitivity study has been conducted for the system with H = 1.0, 
He = 0.304 and ze = 0. Comparing the exponential growth rates 7 obtained 
on various grids, we have determined that the grid with = 60, = 60, 

and No = 96 is sufficient for accurate results. Further increase of the grid 
size by one third in each direction changes 7 by less than 2 %. 

A direct validation of the numerical model is impossible in our case, since 
no accurate and reliable experimental, numerical, or analytical data for the 
system considered in this paper are available. We can, however, rely on 
the extensive verihcation of the numerical method conducted in our earlier 


studies of magnetohydrodynamic shear and convection flows [19|, l2l|, |22, l23|. 
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24 . 25l |. Furthermore, as demonstrated in section 01 the model’s validity 
is supported by the consistency of its predictions with the known data for 
the convection instability in the Rayleigh-Benard layer and in the layer with 
internal heating. 


3. Results 

The prohle of the base-state temperature Tq(z) illustrated in £g. [T]d sug¬ 
gests the possible mechanism of the instability. The temperature field is 
stratihed unstably in the upper part of the electrolyte layer and, to a much 
smaller degree, in the layer A. If sufficiently strong, the stratihcation should 
cause convection flow, which may induce flow in the rest of the battery. In 
this section, we present the results of the computational analysis of this phe¬ 
nomenon. A discussion of the underlying physics and of the implications for 
the battery’s operation are given, respectively, in sections 0] and |5l 

3.1. Detailed analysis of a convection flow 

As a typical example of the instability and the resulting convection flow, 
we consider the case He = 0.304 and Gr = 4.5 x 10®. The evolution of 
the volume-averaged perturbation energies (1231) is shown in Fig. |2l We see 
that after a short initial period the perturbations become dominated by the 
strongest instability mode, and the energies grow exponentially. The growth 
continues till about t = 350, when the levels A' ~ 3 x 10“®, Ee ~ 5 x 
10 “^, corresponding, in our units, to the hnite-amplitude saturation of the 
instability, are reached. After that, the flow’s evolution is nonlinear, with 
slow and irregular oscillations of energies around constant values. 

The perturbation energies computed as in ([23]), but with the averaging 
performed over separated layers A, B, and E, are shown in Fig. [31 We see 
that the amplitude of the perturbations is the highest in the electrolyte and 
the lowest in the bottom layer B. Specihcally, in the developed flow, we hnd, 
for the time-averaged values, Ea/Ee ~ 0.096, Eb/Ee ~ 0.024. 

In order to illustrate the spatial structure of the flow. Figs. 0] and [5] show 
distributions of temperature and velocity in the vertical (drawn through the 
axis) and horizontal cross-sections. 

The exponentially growing dominant mode of the linear instability (see 
Fig. 0|) is a combination of several three-dimensional convection cells located 
mainly within the electrolyte layer. Much weaker motion is generated in 
the top layer A. The motion in the bottom layer B is even weaker. The 
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Figure 2: (a), Volume-averaged kinetic and thermal energies of perturbations (see (IMl) ! 
in the simulation with He = 0.304 and Gr = 4.5 x 10®. (b), Perturbation energies and 
the growth rate 7 (see (l22l) ') during and shortly after the phase of exponential growth. 




Figure 3: Kinetic (a) and thermal (h) energies of perturbations averaged over individual 
layers in the simulation with He = 0.304 and Gr = 4.5 x 10®. 
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convection cells do not form any regular pattern in the horizontal plane and 
have the typical horizontal size about an order of magnitude smaller than 
the diameter of the battery. 

At the nonlinear saturation stage, the convection cells gradually evolve 
into the pattern illustrated in Fig. [5l Individual unsteady cells are still ob¬ 
served, but the flow acquires a degree of axial symmetry. In the electrolyte 
layer, the zones of downward flow are concentrated near the sidewalls and 
around the axis, while the upward flow primarily occurs in a ring at middle 
values of r. The weaker circulation in the bottom and top layers has a similar 
structure, but the opposite circulation sign. 

To evaluate the effect of convection flow and the mixing associated with 
it on the heat transfer in a battery, we compute the two Nusselt numbers 
corresponding to the heat transfer through the top and the bottom walls: 


Nua = 


d{T)/dz\,=H 

dTo/dz\^^^ 


Nub = 


dTo/dz\^^^ 


(47) 


where (T) is the instantaneous temperature field averaged over the entire 
wall surface 

In a steady state, the total heat flux through both walls remains equal 
to the rate of internal heat generation Qq. This implies that the cumulative 
Nusselt number 


Nu 


d{T)/dz\ 

\z=H+ 9 (T)/dz 

'\z=0 

dTQ/dz\ 

\z=H + dTo/dz\^ 

j=0 


1 

2 


{Nua + Nub) 


(48) 


must be equal to one. 

The evolution of Nua and Nub with time for our typical case is shown 
in Fig. [H^. We see that, when a hnite-amplitude convection flow develops, 
the heat transfer through the top wall exceeds the conduction heat transfer 
by about 8%, while the heat transfer through the bottom wall is reduced 
by a similar margin. The cumulative Nusselt number Nu fluctuates around 
unity conhrming that our system’s behavior, while unsteady, corresponds to 
a statistically steady state. 

The asymmetry of the heat flux in the presence of convection can be 
explained by the mixing effect of the convection flow. Fig. |6 ]d shows the 
distribution of the mean temperature T{z), which we obtain by horizontal 
averaging and time averaging over the entire period of the evolution of fully 
developed flow. We see that the change is such that the amplitude of dT/dz 
is slightly increased at the top wall and slightly decreased at the bottom wall. 
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Figure 4: Spatial structure of perturbations during the phase of exponential growth ob¬ 
tained in the simulation with He = 0.304 and Gr = 4.5 x 10^ at t = 300. (a)^ Temperature 
perturbations T — To and velocity vectors (drawn at every second grid point in each di¬ 
rection) in the vertical cross-section ^ = 0 , tt. (b)-(d)^ vertical velocity in the horizontal 
cross-sections z = 0.25 (h)^ z = 0 (cj, and z = 0.75 (d). 
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Figure 5: Spatial structure of developed finite-amplitude convection flow obtained in the 
simulation with He = 0.304 and Gr = 4.5 x 10^ at t = 1500. (a)^ Full temperature T 
and velocity vectors (drawn at every second grid point in each direction) in the vertical 
cross-section ^ = 0, tt. (b)-(d)^ vertical velocity in the horizontal cross-sections z = 0.25 
(b)^ z = 0 (c)^ and z = 0.75 (d). 
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Figure 6: (a), The Nusselt numbers at the top Nua and bottom Nub walls (see (H71) ') and 
the cumulative Nusselt number Nu (see (l48l) ') computed in the simulation with He = 0.304 
and Gr = 4.5 x 10®. (b), Distribution of mean temperature T at the stage of fully developed 
flow for the same parameters. The temperature profile in the base state Tq{z) is shown 
for comparison. 


3.2. Parametric study 

As discussed in section 12.31 we analyze how the convection instability 
changes with the two parameters of the system: the Grashof number Gr 
and the non-dimensional thickness of electrolyte He- The aspect ratio of the 
battery is kept at H = 1, the electrolyte layer is located in the middle of the 
cell, and the Hartmann number is related to Gr by fl2T]) . In every simulation, 
the evolution of the flow is simulated and analyzed in the same manner as for 
the typical example presented in the previous section. We have found that 
in every case, except whose where the base state is stable, the flow evolves 
through the same stages of exponential growth and nonlinear saturation. 

The results are summarized in Table [H which lists all the completed 
simulations and shows the computed quantitative properties: the exponential 
growth rate 7 (see fl2^ ). and the time-averaged perturbations energies E, 
E'^ (see fl23|l h Ea, Eb, Ee, and the Nusselt numbers Nua, Nub (see fH7|) i 
obtained for fully-developed flows. 

The critical Grashof numbers Gr^r, such that the base state is stable at 
Gr < Gr^r and unstable at Gr > Gr^r, can be approximately determined by 
extrapolating the function 'j{Gr) to 7 = 0 and verifying that the base state 
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He 

Gr 

7 

E 


Ea 

Eb 

Ee 

Nua 

Nub 

0.1 

3.0E8 

stable 








0.1 

5.0E8 

0.0005 








0.1 

7.0E8 

0.0022 

2.98E-6 

1.23E-9 

6.49E-6 

3.32E-8 

4.88E-7 

1.02 

0.98 

0.1 

1.0E9 

0.008 

7.26E-6 

4.65E-9 

1.57E-5 

8.28E-8 

1.47E-6 

1.07 

0.93 

0.1 

2.0E9 

0.024 

1.29E-5 

1.29E-8 

2.78E-5 

1.39E-7 

2.80E-6 

1.16 

0.85 

0.1 

3.0E9 

0.033 

1.28E-5 

1.59E-8 

2.76E-5 

1.37E-7 

3.01E-6 

1.19 

0.82 

0.304 

2.0E6 

stable 








0.304 

2.5E6 

0.004 

9.38E-7 

1.27E-7 

2.49E-7 

5.83E-8 

2.73E-6 

1.02 

0.98 

0.304 

3.0E6 

0.013 

1.81E-6 

2.26E-7 

4.68E-7 

1.51E-7 

5.23E-6 

1.04 

0.96 

0.304 

3.5E6 

0.020 

2.30E-6 

3.32E-7 

6.10E-7 

1.82E-7 

6 .66E-6 

1.05 

0.95 

0.304 

4.0E6 

0.027 

2.80E-6 

4.63E-7 

7.63E-7 

2.16E-7 

8.11E-6 

1.06 

0.94 

0.304 

4.5E6 

0.032 

2.85E-6 

4.65E-7 

7.60E-7 

2.17E-7 

8.23E-6 

1.06 

0.94 

0.58 

2.0E5 

stable 








0.58 

2.4E5 

0.0016 








0.58 

2.6E5 

0.007 








0.58 

2.8E5 

0.052 

3.25E-5 

1.76E-5 

2.23E-6 

4.63E-7 

5.78E-5 

1.09 

0.91 

0.58 

3.5E5 

0.07 

4.57E-5 

2.10E-5 

4.89E-6 

6.68E-7 

7.67E-5 

1.12 

0.89 

0.58 

5.0E5 

0.09 

4.75E-5 

3.77E-5 

6.16E-6 

3.46E-7 

7.96E-5 

1.16 

0.84 


Table 1: Summary of the parametric study. The rate of exponential growth 7 (see (l22ll l. 
and the time-averaged perturbations energies E, E'^ (see (l23)l b time-averaged perturbation 
kinetic energies Ea, Eb, Ee, and the Nusselt numbers Nua, Nub (see (1171) 1 of fully- 
developed flows are shown for all the completed simulations. Slight deviations from the 
integral relation between Ea, Eb, Ee, and E and from the relation Nu = 1 (see (1481) 1 are 
due to the time averaging error. In several cases with weak instability (small values of 7 ) 
extension of the simulations beyond the exponential growth stage required very long runs 
and was not conducted. 
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He 

0.1 

0.304 

0.58 

GT (yp 

4.67E8E8 

2.28E6 

2.30E5 

Grl 

4.67E3E3 

5.91E3 

1.51E4 


Table 2: Instability thresholds Gvcr and Grf^. (see (0^) obtained by extrapolation of the 
curves jiGr) from table [T] 

is stable at smaller Gr. The results are presented in Table O We see that, in 
term of our parameters, the batteries become more stable (requiring larger 
Gr for the instability) as the non-dimensional thickness of the electrolyte 
becomes smaller. This can be explained by the fact that the Grashof number 
is dehned using the battery’s radius R as the length scale (see flTTIl h The 
approach is justihed, because it establishes a direct link between the convec¬ 
tion instability and the size and total current of the battery. At the same 
time, it is not consistent with the nature of the instability, which is caused by 
vertical stratification of the base temperature Tq{z) in the electrolyte layer 
and, therefore, has He as the relevant length scale. To address that, the 
table [2] also shows the critical values of the if£;-based Grashof number 

Gr® = Gr/Zi (49) 

We see that these values are closer to each other. An increase of Gr® with 
growing He can be attributed to the increasing constraining effect of the 
sidewalls. This aspect of our results is further discussed in section 01 

The typical structures of the fully developed convection flows in batteries 
with He = 0.58 and He = 0.1 are shown, respectively, in hgures [7] and El 
We see that the principal flow structure remains the same as in the case of 
He = 0.304. It can be described as a pattern of irregular three-dimensional 
comvection cells. Some variations are observed, though. In particular, the 
flow at He = 0.58 is characterized by a smaller number of larger convection 
cells. The flow at = 0.1 demonstrates a flow pattern dominated by sheets 
of strong upward or downward velocity similar to the patterns often found in 
convection in shallow layers (see, e.g. for an example of such a pattern 
in the case of penetrative convection). 

Figures [SHHl and tablelHshow that the partition of the kinetic energy of the 
flow among the three layers changes with He- At He = 0.58, the strongest 
flow is still generated in the electrolyte layer, the flows in the layer A and, 
especially, layer B being much weaker. The situation changes at He = 0.1, 
when the flow in the layer A has the average kinetic energy about an order 
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of magnitude higher than the flow in the layer E. The flow induced in the 
bottom layer B remains weak. 


4. Discussion 


The only plausible explanation of the convection instability reported in 
section [3] is the unstable stratihcation of the base temperature in the upper 
part of the electrolyte layer (see £g. [T]d). This is indicated by the fact that the 
unstable temperature gradient is much higher in the electrolyte layer than in 
the layer A and by the structure of the instability mode illustrated in figure 

SI 

It is interesting to explore the analogy between our case and the convec¬ 
tion instability in a uniformly internally heated layer between two horizontal 
walls (see, e.g. j^) or the classical Rayleigh-Benard instability in a layer 
between stress-free boundaries j^. For that purpose, we use the dimen¬ 
sional height of the unstably stratihed sub-layer H'^/2 as the typical length 
scale and redefine the temperature scale as the typical difference between 
the maximum and minimum base-state temperature in the electrolyte layer 
AT* = Qq /8ke- This results in the effective Rayleigh number dehned 

in terms of our non-dimensional parameters as 


= 


qPePeCe 

ue^-e 


AT *(^1 = 


g(3EpECEQo{H^)^ ^ 


MveI^e 


64 


(50) 


Recalculating from table [H we hnd the instability thresholds Ra^J.^ = 241 
at He = 0.1, Rall^ = 305 at He = 0.304, and Ra^l^ = 785 at He = 0.58. 
The values are of the same order of magnitude as the analogously dehned 
Rall^ = 560 for the internally heated layer Q and Ra^l^ = 277r'^/4 ^ 657 
for the Rayleigh-Benard convection with stress-free boundaries [^ . 

Quantitatively, however, our values are different. Of particular concern is 
the fact that the threshold found for the shallow layer He = 0.1 is about two 


times smaller than in j28| and [2^. This can be explained bt the differences 
in the system’s geometry. The zero shear stress conditions imposed on both 
boundaries in the Rayleigh-Benard case or the no-slip condition on the upper 
boundary in the internal heating case restrict development of convection rolls 
in the unstably stratihed layer and, thus, delay the instability in comparison 
to our case, in which freely movable liquid above and below the unstably 
stratihed layer create ‘soft’ boundary conditions. In order to validate this 
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Figure 7: Spatial structure of developed finite-amplitude convection flow obtained in the 
simulation with He = 0.58 and Gr = 5.0 x 10® at t = 600. (a), Full temperature T and 
velocity vectors (drawn at every second point in each direction) in the vertical cross-section 
9 = 0, TT. (b)-(d)^ vertical velocity in the horizontal cross-sections z = 0.25 (b), z = 0 (c), 
and z = 0.75 (d). 
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Figure 8: Spatial structure of developed finite-amplitude convection flow obtained in the 
simulation with He = 0.1 and Gr = 3.0 x 10® at t = 4488. (a), Full temperature T and 
velocity vectors (drawn at every second point in each direction) in the vertical cross-section 
9 = 0, TT. (b)-(d), vertical velocity in the horizontal cross-sections z = 0.25 (b), z = 0 (c), 
and z = 0.75 (d). 
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Figure 9: Spatial structure of developed finite-amplitude convection flow obtained in the 
simulation with He = 0.304 and Gr = 4.5 x 10® at t = 1500. Vertical velocity Uz and 
velocity vectors (drawn at every grid point) are shown in a part of the vertical cross-section 

0 = 0, TT. 


explanation, we have conducted additional simulations with H = 0.1 and 
He = 0.09, i.e. for a cell, in which the layers A and B are very thin and the 

The 
was 


3 . 
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system approaches the single internally heated layer analyzed in 
instability threshold Ra^J,^ ~ 548, much closer to the results of 
found. 

Flows in the liquid metal layers A and B are induced by the flow de¬ 
veloping in the electrolyte via interfacial coupling. The two major coupling 
mechanisms allowed by our model are the viscous shear stress and conduc¬ 
tion heat transfer across the interfac^. The shear stress is active at both 
the interfaces. It leads to formation of circulation cells in the layers A and 
B that have the same direction of horizontal velocity and, thus, the opposite 
circulation sign as the adjacent convection cell in the electrolyte. In such a 
flow, the zones of upward (downward) velocity in the layer E should approx¬ 
imately correspond to zones of downward (upward) velocity in layers A and 
B. 

The coupling by heat conduction is active only at the A-E interface, where 
upward flows of electrolyte create hot spots, which enhance the unstable 


^There is also electromagnetic coupling, but it is, as we discuss later in this section, is 
weak. 
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stratification in the layer A and may cause convection there. Considered in 
isolation from other coupling mechanisms, this should cause an upward flow 
in the layer A above the zones of upward flow in the layer E. 

Visual inspection of the flow structures in £gures|5]l8]as well as of the addi¬ 
tional illustration in figure [9] suggests that the shear stress is the predominant 
mechanism of interfacial coupling in our system. The heat conduction across 
the interface plays, at best, a secondary role. We note that this conclusion 
as well as the observed general structure of the flow can change when defor¬ 
mation of interface, surface tension, and other so far neglected factors of the 
system’s dynamics are included into the model. 

In accordance with the mechanisms of flow generation, the amplitude of 
convection flow varies among the three layers. It is expected that the primary 
buoyancy-driven flow in the electrolyte layer is stronger than the secondary 
flows in the layers A and B. Table [1] confirms this, but only for the cases 
of thick {He = 0.58) and intermediate {He = 0.304) electrolyte layers. In 
the case of a thin layer {He = 0.1), the flow in the layer A has, on average, 
stronger velocity than the flow in the layer E. This can be explained by the 
restrictive action of close top and bottom boundaries of the layer E, which 
prevents development of large convection cells in this layer. 

The convection introduces mixing that affects the heat transfer through 
the battery. In the framework of our model with insulating sidewalls and top 
and bottom walls kept at the same temperature, the effect is such that the 
heat transfer through the top wall increases, while the heat transfer through 
the bottom wall decreases by approximately the same amount in the time- 
averaged sense (see the values of Nua and Nub in tabled]). The temperature 
peak within the electrolyte layer is suppressed, but not dramatically so (see 
figure (Hb). We conclude that the convection has some, but not very strong 
effect on the temperature regime of the battery. The conclusion may change 
when signihcantly higher Gr are considered. 

Our approach in this study is to take into account the magnetohydrody¬ 
namic (MHD) effects in the form of the Lorentz force Ha^Re~^j x Bq arising 
from the interaction between the electric currents induced by the flow and 
the base magnetic field Bq. This is done to account for possible effects of the 
force on the flow, in particular, for suppression of the unstable convection 
modes, delay of the instability, and formation of anisotropic flow structures 
elongated in the direction of Bq (see, e.g., 0 for an introductory discussion 


3ol | for examples of these effects). The results presented 
so far in this paper, however, do not demonstrate clearly visible MHD el¬ 


and [25|, [20 
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fects. In particular, the anisotropy, which would, in our case, take the form 
of tendency to axially symmetric flow structures, is not seen. 

An explanation to the absence of strong MHD effects is found when we 
calculate the Stuart number 

Y _ _ BlRaE 

Re peU 



which is an estimate of the typical ratio between the Lorentz and inertial 
forces acting in the flow and is typically taken as a measure of the anticipated 
degree of the MHD transformation of the flow. Two corrections are needed to 
hnd the relevant value of N. First, since the MHD effects are much stronger 
in highly conducting metals A and B, conductivity aA = ctr = 10'^cte must 
be used. Second, we see in hgures [5HH] and table [U that the non-dimensional 
velocity of the developed convection flow has small (about 10“^ or smaller) 
amplitude of velocity. This is related to the known fact that the free-fall 
velocity U (see ([5|)) overestimates the actual velocity scale in a convection 
flow. Using U in flSTj) leads to underestimating N and, thus, the strength of 
the MHD effects. 

The effective Stuart number can, therefore, be computed as (we use the 
expressions (ITHil . (|2Tll . (|26il to compute Ha and Re) 


Neff ^ 10^—A = 10®/?2(Yr3/io ^ X (52) 

O-E 

Finally, using the instability threshold values Gver in table |2] we hnd Neff ~ 
4.4 X 10-^ at He = 0.1, 8.9 x lO"® at He = 0.304, and 4.48 x lO"® at 
He = 0.58. We see that Neff <C 1 and, so, the effect of the ma gnet ic held 


on the velocity held of the convection how is negligible (see, e.g., |30|). 


As an approximation, we can take the value Neff = 0.1 as the one, at 
which the MHD ehects become signihcant. According to fl52l) . this corre¬ 
sponds, in the framework of our model, to Gr = 3.35 x 10^®. 


5. Implications for battery design 

We start with the estimates of the critical radius Rcr such that the bat¬ 
teries of larger radii experience the convection instability. Using the values 
of Gver in tableland the formula fl27|l we hnd Rcr = 3.6 cm at He = 0.1, 
Rcr = 1-2 cm at He = 0.304, and Rcr = 0.78 cm at He = 0.58. 
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We can also estimate the expected intensity of the convection flow. Par- 


ticnlarly interesting is the flow in the bottom layer. As discnssed in [4, [15 


this flow plays a potentially benehcial role, as it enhances the mixing of metal 
B and componnd A-B, thns redncing a limiting factor of the reaction rate. 
The typical valne of the non-dimensional velocity can be evalnated from ta¬ 
ble [Has Ub ~ . Mnltiplying by the velocity scale U = see 

dSD and nsing the physical properties from section 12.31 and cnrrent Jq = 1 
A/cm^, we hnd the typical dimensional amplitnde of velocity in the bottom 
layer as 

~ m/s. (53) 


This gives rather low valnes of velocity. For example, in the particnlarly 
interesting case of thin electrolyte layer = 0.1 at Gr = 3 x 10®, we 
hnd f/f ™ ~ 0.3 mm/s. Even lower valnes are considered in the other cases 
considered in onr stndy. We mnst note, however, that the example above 
corresponds to a small cell (radins abont 5.2 cm with the cnrrent and physical 
parameters nsed in onr work). Mnch stronger circulation and turbulent how 
are expected in larger batteries. For small batteries, desired mixing can 
be produced by other mechanisms, such as bottom heating, electrovortex 
instability, etc. 

We can make preliminary estimate of the interface deformation cause by 
the convection how. The question is interesting for two reasons. First, the 
possibility of signihcant deformation would invalidate our model, based on 
the assumption of non-deformable interfaces. Second, a deformation so string 
as to cause rupture of the electrolyte layer at some point is highly undesirable, 
since this would lead to short circuit between the metal layers and, thus, a 
disruption of the battery’s operation. The effect cannot be analyzed directly 
within the framework of our model, but can be approached via a qualitative 
energy arguments similar to those used for the Tayler instability in [l^ . We 
recall that the densities of the melts in a real battery are strongly different and 
consider the practically most interesting case of thin electrolyte layer. The 
estimate is based on comparison between the kinetic energy of the how with 
the gravitational potential of the interface deformation. Since the density pb 
of the melt in the bottom layer B is much higher than pe and pa and since 
the how is the strongest in the layer A, we only consider the deformation of 
the upper interface. A displacement of a huid particle at the interface by the 
dimensional distance h* increases the specihc gravitational potential by 


Epot = {pE - Pa) gh*, 


(54) 
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which should be of the same order of magnitude as the specific kinetic energy 
of the flow 




(55) 


where U* is the dimensional mean velocity in the layer A. This leads to the 
estimate 

1 PA !/•" 

tl ~- 

2 pe - pA g 

or, in the non-dimensional form 



h = 



1 Pa 

2 Pe — Pa 


Fr\ 


(57) 


where H'^ is the dimensional thickness of the undisturbed electrolyte layer 
and 


Fr 


U* 


(58) 


is the Fronde number. 

The event of short circuit between the metal electrodes corresponds to 
h* > i.e. the Fronde number exceeding the critical value 


Ft 

± / cr 




(59) 


The typical values of the coefficient {pE — Pa) / Pa is about 2. We consider 
the case of He = 0.1 and Gr = 3 x 10®, where the strongest fiow is found 
above the thinnest electrolyte layer, and estimate the square of the mean 
velocity as U*‘^ ~ U^Ea, where U is the free-fall velocity (jS]) and Ea = 
2.76 X 10“® is the non-dimensional kinetic energy density from table [T] This 


gives 

fa ~ EaFL = 

g^E CTeKe 


(60) 


With the physical parameters used in our study, Gr = 3 x 10® corresponds 
to the battery radius R = 5.2 cm, which gives ~ 3 X 10 We see 
that the amplitude of the deformation of the interface is about four orders of 
magnitude smaller than the thickness of the electrolyte layer. The validity 
of our modeling assumption is confirmed, and the danger of the rupture 
of interface by the convection ffow should be seen as immaterial. These 
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conclusions are, of course, only valid for the range of parameters considered 
in our paper. The situation can be completely different in larger battery 

cells. 

Finally, we estimate the size of the battery, at which the convection flow 
is signihcantly affected by the magnetic held Bq. Following the calculations 
at the end of the previous section we hnd, for Gr = 3.35 x 10^®, the radius 
R = 1.33 m. At smaller sizes, the effect can be ignored in the analysis. At 
the same time, attention may be needed to the effects of the magnetic helds 
generated externally, for example, by the current supply lines. 

6. Concluding remarks 

The main conclusion of the analysis presented in this paper is that thermal 
convection is virtually inevitable in liquid metal batteries. It is present in 
even small-scale prototypes. In the practically most interesting case of thin 
electrolyte layer, the how is the strongest in the top layer containing the 
liquid metal and the weakest in the bottom layer. 

The ehect of the how on the battery’s operation is found to be weak. The 
velocity in the bottom metal layer is too small to generate strong mixing of 
reactants. No signihcant deformation of the electrolyte-metal interfaces is 
expected. 

These conclusions are, however, only valid for small (a few cm in radius) 
batteries, for which computations have been conducted in this study. The 
ehect can and should be anticipated to be mush stronger in larger batteries, 
for which the Grashof number Gr ~ is many orders of magnitude larger. 

It appears that the most interesting direction of future work is the anal¬ 
ysis of operation of large {R ~ 0.1-1 m) batteries. In addition to being 
computationally challenging, this will require more complex physical mod¬ 
els, which will include deformation of interfaces, dynamics of the base current 
distribution Jq, etc. Ultimately, it will be necessary to develop a comprehen¬ 
sive model that describes the multiple mechanisms of instability and how 
generation as coupled physical phenomena. 
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